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that, with quadrature methods each determinant can be computed separately using the operators 
£L a = D\D a and D./, = D' b Df,. We also discuss using bootstrap re-sampling to remove the bias 
from the determinant estimator. 



XXIX International Symposium on Lattice Field Theory 
July 10-16, 2011 

Squaw Valley, Lake Tahoe, California 



* Speaker. 

' Current address: Computation-based Science and Technology Research Center (CaSToRC), The Cyprus Institute, 
15 Kypranoros Str., RO. Box 27456, 1645 Nicosia, Cyprus 



© Copyright owned by the author(s) under the terms of the Creative Commons Attribution-NonCommercial-ShareAlike Licence. 



http: //pos . sissa . it/ 



Application of Quadrature Methods for Re-Weighting in Lattice QCD 



Abdou Abdel-Rehim 



1. Introduction 

Re-weighting is a useful tool that has been applied recently in lattice QCD simulations in 
different contexts. The main idea is that the expectation value of an observable O with respect to an 
action S b can be written in terms of expectation values with respect to another action S a as follows: 



(0) b := ^ J @ U W e ~ 

Z JL± [ @U 0[U)w[U]e- 

(1.1) 



Z b . 
Za 

(Ow) a 
{w)a 



Si, 



where Z a and Z b are the corresponding partition functions, U is the gauge configuration and w[U] is 
called the reweighting factor corresponding to the gauge configuration U and is defined by w[U] := 
e l^u] ■ The labels "a" and "b" refers to a particular choice of the action parameters such as the sea 
quark mass. Equation [HI] is exact. However, when using importance sampling and for a finite set 
of gauge configurations, reweighting could be biased. Because of this bias, reweighting is only 
reliable when there is a small change between the actions S a and S b , and the reweighting factors 
normalized by the ensemble mean are 0(1). 

An important application of reweighting methods is mass reweighting in which observables 
corresponding to a sea quark mass m b are computed using configurations generated with a sea 
quark mass m a . In this case the reweighting factor is given by 

det[(D(U;m b yD(U;m b ))^} 
w[U\ — w— , (i.z) 

det[(D(U;m a yD(U;m a ))-} 

where D(U;m) is the Dirac operator for gauge configuration U and bare quark mass m and tif is the 
number of degenerate flavors being re-weighted. Applications of mass reweighting include tuning 
the strange quark mass in dynamical Domain-Wall simulations to its physical value [gj, computing 
the strange Nucleon sigma term using the Feynman-Hellman theorem [Q], simulating the e regime 
[^], and tuning to the physical point in 2+1 simulations of the PACS-CS collaboration 

A common approach for computing the ratio of determinants is using pseudofermions. This 
has been the main technique used for mass reweighting [||, |[ [5[]. In this approach, the relation 

1 f^^e-5 tQ « 

- ' (1.3) 
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is used to compute the determinant of a matrix O. The right-hand side of Equation ( 1.3) is computed 



stochastically as an average over random gaussian fields ^ generated with the distribution e % f( = as 
where the notation < g >^ means expectation value of g over For a reliable estimate of det(£2) 



using Equation ([L4J), £2 need to be close to the identity matrix. For nf degenerate flavors, 



a=[D a (DlD h )- l Dlfl. (1.5) 
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For small changes in the fermion action, Q. will be close to the identity matrix. In addition, it 
is possible to improve the estimation of det(H) by dividing the change of the Dirac operator into 
smaller steps and writing the ratio of determinants as a product of k + 1 factors 

det(A,) = det(D fl ) det(Di) det(D fc _i) det(D t ) 
det(Di) det(Di) det(D 2 ) "' det(D t ) det(D,,)' 

such that the ratio of determinants in each factor is close to one. The pseudofermion approach 
has the advantage of giving directly an unbiased estimator of the ratio of determinants. It however 
requires the solution of a linear system for each Another difficulty with the pseudofermion 
approach is that, for a single flavor reweighting, it requires using a rational approximation for the 
square root of a matrix which in turn requires a multi-shift solver. We study here an alternative 
way of computing the ratio of determinants based on the relation det(£2) = e trace [ lo g(^)]_ i n 0U r 
approach, the trace of log(O) is evaluated using noise methods and the matrix elements of log(O) 
between noise vectors are computed using Gauss quadratures This approach has been proposed 
in the context of generating dynamical electromagnetic field configurations on top of existing QCD 
configurations [ph. In this paper, we elaborate on using this approach also for mass reweighting. 
In addition, we show that the need to solve a linear system for each noise vector can be avoided 
by computing the determinant of each Dirac operator separately and then take the ratio. This leads 
to a faster evaluation of the reweighting factor than what is usually used in the literature. Another 
advantage of using the relation det(£2) = e trace [ 1 °g(^ 1 )] j s that including a fractional power is trivial 
since 

det(r a ) = £ trace [ lo g( r ")] = e «trace[log(r)]_ q 

Finally, variance reduction techniques such as breaking the ratio of determinants into factors and 
dilution can be used. It is noted that dilution couldn't be used with pseudofermions. 



2. Quadrature and Noise approach 

Let z be a vector whose elements are random variables satisfying 

<ZiZj>=<z*z j >=8 ij . (2.1) 
Such vectors are called noise vectors. For a matrix H, we have 

(z i Hz)noise = trace[H], 

variz 1 HzUise = £ \Hij\ 2 + H i} H* fi , (2.2) 

¥J 

where (...)noise and var(. . .) no i se means the expectation value and variance over the set of noise 
vectors. For a positive definite Hermitian matrix H, the element vHz can be written as a Riemann- 
Stieltjes integral which can then be computed using Gauss quadrature methods (see [^] for more 
details). For computing the ratio of determinants, one could define Q. = D a (DtDb)~ l Dl and take 



H := log(Q.), then use Equation [1.7]. In this case, Q. is close to the identity, however, each 
application of II to a vector requires the solution of a linear system which is expensive. We call 
this the standard method. Alternatively, one could compute each determinant separately using the 
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operators £l a = D\D a and £2/, = D b D\,. Although these operators are not close to the identity, they 
don't involve the inverse making this approach potentially faster. We call this the difference method. 

Our tests are done using Clover fermions. In this case it is advantageous to use use even-odd 
preconditioning. The Dirac operator is written in the form 

D= (D ee D eo \( 1 0\ f G 0\ f 1 D~ l D eo \ 
\D oe D 00 j \D oe D-} I J \ QJ \0 1 )' 

where G = D ee and Q = D u0 — D oe D~ e l D eo is the Schur complement. The determinant of the 
Dirac operator is then given by det(D) = det(G)det(<2). For this decomposition, det(G) can be 
computed exactly and cheaply while det (Q) is estimated using noisy estimators. 

Noisy estimators give an unbiased estimate for the trace of the logarithm. In order to get 
an unbiased estimate of the determinant by taking the exponential of the trace of the logarithm, 
additional techniques are needed. There are two possible ways to obtain such unbiased estimator. 
The first method is using a resampling technique such as bootstrap (or jackknife). The second 
method, which is applicable to functions given as a power series, is based on stochastic summation 
of the series In our tests, we used bootstrap to obtain an unbiased estimate. Let^i ,x%, . . . ,Xn be 
a sample of size N of measurements of a random variable x giving a sample average x = 4 Y!i=\ x i- 
To obtain an unbiased estimator of the function g((x)) where (x) is the true mean we generate Nb 
bootstrap samples x} from the original data where i = 1,2, ... ,N and r\ = 1,2, . . . ,Nb- Compute 

1 N i n b i N B 

The bootstrap estimator for g({x)) is unbiased up to ^(ttt) and its error is given by 



*«*»«2 g (50-*» a^ rap '= y / ^ T V V) 2 : - / • (2-5) 
The bootstrap approach is applicable to a general function g. In our case g(x) = e x . 



3. Results 

The method is tested on 16 3 x 48 lattices with tadpole improved tree-level Symanzik gauge 
action at /3 = 6.5 corresponding to lattice spacing a 0.09 fm. A tadpole improved clover fermion 
action is used with one level of stout smeared links with smearing parameter p = 0.125. We have 
three degenerate flavors with quark masses roughly in the range of the strange quark mass. Con- 
figurations were generated with bare quark mass parameters m q = —0.170,-0.175,-0.180 with 
corresponding lattice pion masses am n = 0.3570(16), 0.3302(25), 0.3097(18) respectively. The 
corresponding pion masses in physical units are m 7l (MeV) = 780(35), 720(55), 680(39). In our 
tests, reweighting was used to compute observables for bare sea quark masses m q = —0.175, and 
—0.180 by reweighting configurations generated with m q = —0.170. We have analyzed 700 con- 
figurations. The fact that we have configurations generated with m q = —0.175 and —0.180 allowed 
us to check the conectness of the reweighting procedure by comparing to expectation values of 
observables computed on configurations generated with the correct sea quark mass. Calculations 
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done using Chroma software [Qj. In order to improve the accuracy of the computation of the matrix 
element z*log(Cl)z, we use double precision for the quadrature Lanczos algorithm with single pre- 
cision for the matrix-vector multiplication. We use complex Z4 noise vectors where each element 
of z can have the values 1,-1,/, — i with equal probability. 

First, we compare trace (log (<2j ( <2 a )) — trace (log (Q^ Qb)) (difference method) to 
tmce(log{Q a (QlQb)~ l Ql\) (standard method). Mathematically, the results should be identical and 
stochastically they should be consistent within errors. For this test, we use 1000 Z4 noise vectors 
and m a = —0.170 and nib = —0.180. To investigate the effect of numerical precision, we do the 
comparison for the situation where all the computation is done in double precision (both matrix- 
vector multiplication, the linear solver, and the Lanczos quadrature algorithm is done in double 
precision), and when mixed precision calculation is used (matrix-vector multiplication and the 
linear solver is done in single precision, while dot products inside the Lanczos quadrature algorithm 
are done in double precision). The Conjugate Gradient (CG) algorithm is used to solve linear 
systems when the standard method is used. In the case of double precision calculation, the tolerance 
for the linear system is set to le — 10 and the quadrature for computing the matrix element is also 
computed to converge to the same tolerance. In the mixed precision calculation, the linear solver 
and the quadrature is set to converge to a tolerance of le — 7 when the standard method is applied. 
When the difference method is applied, the quadrature for each matrix element is set to le — 9. The 
reason it was necessary to set the tolerance of the quadrature in the case of the difference method 
to le — 9 instead of le — 7 as it was in the standard case is to have the same number of significant 
digits (7 digits) after taking the difference. In Table [I], we compare the results which show that 
we get results which are consistent with each other within statistical errors. This is important 
since the difference method is much faster than the standard method. In Table ||, we compare the 
cost of each method in terms of the number of matrix-vector multiplication each method uses. In 
the case of the standard method, the Lanczos quadrature algorithm takes few iterations, about 5 
iterations, however each iteration involves the solution of a linear system using CG which takes 
about 500 iterations (a matrix-vector product here means multiplication of a vector with Q?Q). 
The quadrature algorithm converges in few iterations because the compound operator is close to 
the identity. On the other hand, for the difference method, no linear system need to be solved, 
however, the Lanczos quadrature algorithm takes more iterations to converge since the operators 
QlQa and QlQb are far from the identity. The total cost of the difference method is much smaller 
than the standard method. 

Since both methods give the same results, we use the difference method for mass re-weighting. 
We reweight configurations generated with sea quark mass —0.170 to compute observables corre- 
sponding to sea quark masses —0.175 and —0.180. We computed the reweighting factor for 700 
configurations and for the production calculations we use only 500 noises and a mixed precision 
calculation in which the matrix-vector products are done in single precision while dot products in 
the Lanczos quadrature algorithm are done in double precision. The matrix elements were made to 
converge to tolerance le — 9. 

In Figure [TJ, we show the normalized re-weighting factors. As shown in the figure, there is 
a larger fraction of the configurations with small reweighting factor for the —0.180 case than the 
—0.175 case. Having the re- weighting factor available, we can now re-weight any physical ob- 
servable of interest. As an example, we look at re-weighting the average thin plaquette. In Figure 
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|2[ we compare re-weighted various thin plaquettes to the correct values obtained from configura- 
tions generated with the correct sea quark mass. The results show that 700 configurations provide 
enough statistics to re-weight from sea quark mass —0.170 to —0.175, but not to correctly reweight 
from -0.170 to -0.180. 



Config. No. 


Standard Method 
Double Precision 


Difference Method 
Double Precision 


Standard Method 
Mixed Precision 


Difference Method 
Mixed Precision 


1 


5485.76(9) 


5485.78(8) 


5485.91(9) 


5485.77(8) 


2 


5484.34(9) 


5484.35(8) 


5484.49(9) 


5484.33(8) 


3 


5486.57(9) 


5486.63(8) 


5486.73(9) 


5486.62(8) 



Table 1: Comparison of trace{log[g a (g^gf,) _1 gJ]} (standard method) and trace{log[gJg„]} — 
trace{log[£>i2f>]} (difference method) on three configurations where Q a corresponds to m a — —0.170 and 
Qb corresponds to = —0.180 using 1000 noises without dilution. 



Config. Num. 


Standard Method 
Mixed Precision 


Difference Method 
Mixed Precision 


Standard Method 
Double Precision 


Difference Method 
Double Precision 


1 


5+2671 


105+123 


6+5164 


137+162 


2 


5+2652 


103+121 


6+5007 


137+163 


3 


5+2685 


106+123 


6+5172 


138+165 



Table 2: Comparison of the number of matrix-vector products used by the standard and difference methods 
for a single noise vector when S a corresponds to mo = —0.170 and corresponds to mo = —0.180. For the 
standard method, the first number is the number of Lanczos algorithm iterations and the second number is 
the total number of iterations used by the CG solver during these iterations. For the difference method, the 
two numbers correspond to the number of iterations used by each Lanczos algorithm. 




Configuration Number Configuration Number 



Figure 1: The normalized re-weighting factors. Left: reweighting from sea quark mass —0.170 to —0.175. 
Right: reweighting from sea quark mass —0.170 to —0.180. 



4. Conclusions 

The application of Lanczos quadrature algorithm in mass re-weighting is studied. It is shown 
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average thin plaquette, m =-0.170 — > -0.175 



average thin plaquette, m =-0.170 — > -0.180 



m sea =-0.175 
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Figure 2: Reweighted thin plaquettes compared with correct values. The horizontal axis corresponds to 
different types of thin plaquettes (different planes, plaquettes in the spatial directions only, etc.). Left: 
reweighting from sea quark mass —0.170 to —0.175. Right: reweighting from sea quark mass —0.170 to 
-0.180. 

that it is possible to avoid the need to solve a linear system and apply the algorithm to compute each 
determinant separately before taking the ratio. We also showed how to get an un-biased estimator 
for the reweighting factors using bootstrap method. Having configurations generated with the sea 
quark mass we are trying to re-weight to, allowed us to compare the re-weighted observable to the 
correct ones and check for the range of applicability of the re-weighting procedure. 
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